A new approach to the integration of rotational 
motion in systems with interacting rigid bodies 

Igor P. Omelyan 

Institute for Condensed Matter Physics, National Ukrainian Academy of Sciences, 
1 Svientsitsky st., UA-290011 Lviv, Ukraine. E-mail: nep@icmp.lviv.ua 



(February 2, 2008) 



MSG numbers: 65C20; 65-04; 68U20; 70E15; 82A71 

Keywords: Numerical methods; Motion of rigid bodies; Long-time integration 
Molecular dynamics simulation; Polyatomic molecules; Machine computation 



1 



A running head: 
Numerical integration of rigid-body motion 



Dr. I. P. Omelyan 

Institute for Condensed Matter Physics 
National Ukrainian Academy of Sciences 
1 Svientsitsky st., UA-290011 Lviv 
UKRAINE 

Tel/Fax: +380-322-761978 
E-mail: nep@icmp.lviv.ua 



2 



Abstract 



A new approach is developed to integrate numerically the equations 
of motion for systems of interacting rigid polyatomic molecules. With the 
aid of a leapfrog framework, we directly involve principal angular veloc- 
ities into the integration, whereas orientational positions are expressed 
in terms of cither principal axes or quaternions. As a result, the rigid- 
ness of molecules appears to be an integral of motion, despite the atom 
trajectories are evaluated approximately. The algorithm derived is free 
of any iterative procedures and it allows to perform both energy- and 
temperature-conserving simulations. The corresponding integrators are 
time reversible but the symplectic behavior is only achieved in mean. 
Symplectic versions are also described. They provide the conservation of 
volume in phase space precisely at each time step and, moreover, lead 
to exact solutions for angular velocities in the inertial-motion regime. It 
is shown that the algorithm exhibits excellent stability properties and 
conserves the energy even somewhat better than the atomic-constraint 
technique. 
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I. INTRODUCTION 



A lot of theories in physical chemistry treats the hquid as a system of classical rigid 
bodies. The method of molecular dynamics (MD) remains the main tool for investi- 
gation of such a model. An important problem in MD simulations is the development 
of stable and efficient algorithms for integrating the equations of motion with orienta- 
tional degrees of freedom. The straightforward parameterization of these degrees, Euler 
angles, is very inefficient for numerical calculations because of singularities inherent in 
the description [1-3] . Within singularity free approaches, the orientations of molecules 
typically are expressed in terms of quaternions [4-6]. High-order Gear methods were 
utilized to integrate the rigid-body equations of motion in early investigations [7-10]. 
These schemes are not reversible or symplectic, and it is not clear that the extra or- 
der obtained is relevant, since they quickly become exhibit poor long-term stability of 
energy with increasing the step size [5, 11]. 

Various alternatives to the Gear approach have been proposed and implemented 
over the years. These include Verlet [12], velocity Verlet [13], leapfrog [14], and Bee- 
man [15] integrators, which are based on the Stormer central difference approximation 
[7, 16] of accelerations. Such integrators proved to be the most efficient, because high 
accuracy can be reached with minimal cost measured in terms of force evaluations. 
The main problem with these methods is that they were initially constructed, in fact, 
to integrate translational motion assuming the velocity-independence of the accelera- 
tions, and, therefore, additional modifications are necessary to use them for rotational 
dynamics. 

In the atomic approach [17-21], the parameterization of orientational degrees of 
freedom is circumvented by involving individual Cartesian coordinates of atomic sites. 
As a result, the Verlet algorithm and its velocity versions appear to be directly appli- 
cable. The dynamics is determined by integrating the equations of motion for these 
sites, subject to constraints that the intramolecular bond distances are fixed. Un- 
til now, the great majority of MD simulations on polyatomics are performed using 
the atomic-constraint technique due to its exceptional long-time stability. Although 
this approach can be applied, in principle, to arbitrary systems including the case of 
flexible molecules, it has some disadvantages. For example, to exactly reproduce the 
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rigid molecular structure, complicated systems of nonlinear equations are needed to 
be solved by iterations at each time step of the integration process. In general, the 
convergence is not guaranteed [22] and looping becomes possible already at relatively 
small step sizes, especially for macromolecules with bond angle reactions. An addi- 
tional complexity arises for point molecules with embedded multipoles, since then the 
forces are not easily decomposable into direct site-site contributions. 

In order to obviate these difficulties, Ahlrichs and Erode have devised a hybrid 
method [23] in which the principal axes of molecules are considered as pseudo particles 
and constraint forces are introduced to maintain their orthonormality. The principal 
axes were evaluated within the Verlet framework via a recursive procedure which does 
not solve exactly the constraint equations to convergence, but instead writes the rota- 
tional matrix as an exponential for sum of some anti-symmetric matrices, restricting 
by a finite number of terms. It was established, however, that the exponential prop- 
agation leads to much worse results on energy conservation than those obtained in 
the atomic-constraint technique. Recently, acting in the spirit of the pseudo-particle 
formalism, a leapfrog-hke algorithm has also been proposed [24]. In this algorithm 
the entire rotation matrix and the corresponding conjugate momentum are treated as 
dynamical variables, and the matrix of constrain forces is evaluated exactly. In such a 
way, the lost precision was reproduced, but this required again to find iterative solu- 
tions to highly nonlinear equations. Moreover, since velocities do not appear explicitly, 
it is hard to extend the pseudo-particle approach to thermostatted versions. 

The first efforts on adopting the Stormer group of integrators to rotational motion 
in its pure form were done by Fincham [25-26]. As a result, explicit and implicit 
angular-momentum leapfrog algorithms have been introduced. In the case of a more 
accurate implicit version, the system of four nonlinear equations per molecule is solved 
iteratively for the same number of quaternion components [26]. As was soon realized 
[27], the Fincham's algorithms are not very efficient in energy-conserving simulations, 
given that the artificial rescaling procedure is used to maintain the unit quaternion 
norm, and additional transformations with approximately computed rotational matri- 
ces and angular momenta are necessary to evaluate the principal components of angular 
velocities. The question of how to replace the crude renormalization by a more natural 
procedure has been considered too [28]. As a consequence, the quaternion dynamics 
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with constraints was formulated within the velocity Verlet framework. Similar ideas 
were used to adopt the pseudo-particle approach to velocity- Verlet and Beeman integra- 
tors [27] . It has been shown that such algorithms conserve the total energy better with 
respect to the implicit leapfrog integrator [26] but worse than the atomic-constraint 
method, especially in the case of long-duration simulations with large step sizes. 

Among other techniques developed in recent years it is necessary to mention the 
multiple time scale (MTS) integration [29-34]. The reversible reference system prop- 
agator (RESP) algorithm by Berne at al. [32] is a general approach which yields a 
whole family of MTS integrators. The basic idea of this approach lies in using a short 
time step for the rapidly varying fast forces, and a large time step for the nonbonding 
slow forces arising from the interactions at long interparticle distances. Then the most 
time-consuming slow forces will be recalculated less frequently than in usual methods, 
saving significantly the computer time. It is worth remarking that the RESP approach, 
being originally constructed on Cartesian coordinates, can be adapted [35] to orienta- 
tional variables. There are some arguments to stay that a reformulation of the RESP 
approach within these variables can be more convenient for applications (see Sec. IV) . 
Prom this point of view, and taking into account that the fastest motions are handled 
(inside the innermost cycle with the least step size) using standard algorithms, the de- 
veloping of methods for direct integration of orientational variables presents an interest 
in the context of the MTS approach as well. 

Quite recently, to improve the efficiency of integration in orientational space, a new 
angular- velocity leapfrog algorithm has been proposed [36]. Its main advantages were 
the intrinsic conservation of rigid structures and high stability properties. However, a 
common drawback, inherent in all precise long-term integrators on rigid polyatomics, 
still remained, namely, the necessity to solve by iteration nonlinear equations. 

In the present paper we develop the angular-velocity algorithm by avoiding any 
iterative procedures. A thermostatted version is also introduced. We demonstrate 
that the integrators derived are time reversible, symplcctic and the total energy is 
conserved even better than within the atomic-constraint method. It is discussed that 
our approach can be especially useful for simulations of systems at high temperatures, 
where the inertial-motion regime plays an important role in rotational dynamics. 
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II. ADVANCED LEAPFROG ALGORITHM 



We shall deal with a classical collection of N rigid molecules composed of M point 
interaction sites. According to the molecular approach, the dynamics for such a sys- 
tem can be considered in view of translational and rotational motions. The trans- 
lational motion is expressed in terms of the center of mass r j = Y,^=i ^a^t /m oi 
each molecule (i — 1, . . . ,N), where denotes the positions of atomic site a within 
molecule i and m — f^a and nia are the masses of a separate molecule and 

partial atoms, respectively. Then the translational part of time evolution can be deter- 
mined by the first-order differential equations: drj/di = Vj and mdvj/dt = fj, where 
fj = Yl!j-'ah ~'^^\) the force acting on molecule i due to the site-site interactions 

if^ with the other {j ^ i) molecules and Vj designates the center-of-mass velocity. 

A. Rotational motion in body-vector and quaternion representations 



In the body- vector representation [4, 23, 27], the Cartesian coordinates of three 
principal axes XYZ of the molecule are assumed to be orientational variables. These 
variables can be collected into the 3x3 orthonormal matrices Aj, so that the site posi- 
tions in the laboratory frame are: r^{t) — ri{t)+Af{t)Aa, where Aa = (A^, Ay, A^)"^ 
is a vector-column of these positions in the body frame attached to the molecule 
and A+ denotes the matrix transposed to A. The time- independent quantities Aa 
(a = 1, . . . , M) are completely defined by the rigid molecular structure. The rate of 
change in time of orientational matrices can be given as 



dA, 
dt 
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where Jl^, and Jl^ are principal components of the angular velocity fij, and W is 
a skewsymmetric matrix associated with fij, i.e., W"^(fij) = — W(f2i). 
In an alternative approach [4, 9], the matrix Aj is a function. 
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of the four-component quaternion = {Ci, Vi^ Xi)'^ ■ The time derivatives of can 
be presented [36] in the form 
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(3) 



where Q(r2j) is a skewsymmetric matrix again. The normahzation condition q^^ = 
+ Vi + Ci + Xi — which ensures the orthonormahty AjA^ = I of Aj = A(qj), 
where I designates the unit matrix, has been used to obtain Eq. (3). 

The relations (1) and (3) for orientational coordinates need to be supplemented by 
the Euler's equations for angular velocities 

dfiV 



dt 
dQi 



K'x + {Jy - Jz) flV^ 



Y- 



dt 

dfi*2 

'~dt 



K'y + {Jz-Jx)^'z^\. 



^K'z + {Jx-Jy)^'x^V, 



(4) 



where Jq denote the time-independent moments inertia of the molecule along its prin- 
cipal axes {a — X, Y, Z) and K]^ are body-frame components, Kj = Ajkj, of the torque 
= Y^'j-^^i^l ~ '^d^^ij exerted on molecule i with respect to its center of mass. 



B. Evaluation of angular velocities and coordinates 

Let {rj(t), Vj(t — |), Si(t), f2j(t — |)} be an initial spatially-velocity configuration of 
the system, where the velocities and positions are defined on alternate half-time steps 
with h being the fixed step size and Sj(i) = Aj(t) or qj(i) are the orientational coor- 
dinates for the principal- axis or quaternion representations, respectively. The transla- 
tional variables can be integrated applying the usual [14] leapfrog algorithm 



v,(i + I) = Yi{t - I) + Mi{t)/m + 0{h^) , 
Yi{t + h)= ri{t) + hvi{t + I) + 0{h^) 



(5) 



in which forces fi{t) are explicitly evaluated in terms of known spatial coordinates rj(t) 
and Si{t). As can be verified easily, expanding the left- and right-hand sides of both the 



lines of Eq. (5) into Taylor series over h, the algorithm produces truncation single-step 
errors of order in coordinates and velocities. 

In the case of rotational motion it is not a simple matter to adopt the standard 
leapfrog scheme, since the principal angular accelerations are velocity-dependent and 
the time derivatives of orientational coordinates depend not only on the angular velocity 
but also on the coordinates themselves. These problems were handled previously by 
Fincham [26] in his angular-momentum versions of the leapfrog algorithm. Recently 
[27, 36], it was shown that more superior techniques follow when principal angular 
velocities are involved directly into the integration. Within the leapfrog framework, 
mid-step values for the angular velocities can be evaluated by writing 

Kit + 1) = n^t - 1) + h[Kit) + (j^ - J,) ^rf,m\it)]/Ja + o{h') , (6) 

where Euler equations (4) have been taken into account, 7) denote an array of 

three cyclic permutations for (X, F, Z), and torques Kl^it) are computed via coordinates 
Yi{t) and Si(i). Equation (6) presents a rotational-motion analog for the first hue of (5) 
but it must be complemented by an interpolation of the products of angular velocities 
to on-step levels of time. It is quite naturally to perform such an interpolation in the 
form 

m)^.^) = I Im - - 1) + + iHit + D] + ^(^') ' (7) 

where 0{h?) uncertainties are in the self-consistency with the second-order accuracy 
of angular- velocity propagation (6). In view of (7), vector expression (6) is an implicit 
equation with respect to r2j(t + |), which allows to be solved by iteration [36]. 

Similarly to the translational-position evaluations (second line of Eq. (5)), we inte- 
grate orientational coordinates, 

Si{t + h)^ Si{t) + hHi{t + |)S,(t + I) + 0{h^) , (8) 

for the cases of principal- axis vectors (Sj = Aj, Hj = W^) and quaternion (Sj = 
qj,Hj = |Qj) representations, where Eqs. (1) and (3) have been used. The matrices 
Wj = W(f2i) and Qj = Q(f2i) are calculated in (8) using the already defined values 
of r2i(t + |), whereas the obvious choice for mid-step values of orientational variables 
is 
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s^(^ + 1) = i [Ht) + s^(^ + h)] + 0{h^) . (9) 

Equations (8) and (9) constitute, in fact, a system of linear equations with respect to 
elements of Aj(t + h) or qj(t + h), which, is solved analytically, 

Si{t + h)^[I- |H,(i + + |H,(t + |)]S,(t) + 0{h') , (10) 

where it is understood that I designates either three- or four-dimensional unit matrix 
in the principal- axis or quaternion domains, respectively. 

Taking into account expressions (1) and (3) for matrix Hj, the result (10) can be 
written more explicitly, 

*(* + = ''^"flol'j'jt''^' - ft) q.(«) , 

16 « V* 2/ 

where [Pijo/? = i^of^^ denotes a symmetric matrix which, as the matrices W, and Qj, 
is computed in (11) using the angular velocities at middle- step time t + |. It can be 
checked that the matrix (I — eH)~^(I -|- eH) is orthonormal at arbitrary values of e, 
provided H+ = — H. As a result, the 3x3 and 4x4 evolution matrices Dj and Gj 
are orthonormal by construction as well. Using the equahties W| = — Qfl and 
= — ^^I, these matrices can be presented in the exclusive compact form 

Dj(i, h) = exp[(^jWi/04+| > ^ii^^ M = exp[0jQj/Oj]^_^| , (12) 

where cpi = arcsin[/ia(t + |)/(l + f + 'Pi = arcsin[/i(]i(t + |)/(2 + f l]2(t+|))] 
and exp designates the matrix exponential. Then it becomes clear that matrices Dj and 
Gi define three- and four-dimensional rotations on angles ipi and 4>i in the laboratory 
frame and quaternion space, respectively. In the first case the rotation is performed 
around the unit vector fij/Qj|j^h,, whereas in the second one it is carried out around a 
virtual orth which is perpendicular to all four orths of the quaternion space. 

From the above, the following important statement emerges immediately. If initially 
the orthonormality of Aj and the unit norm of are fulfilled, they will be satisfied 
perfectly in future, despite an approximate character of the evaluated trajectories. This 
excellent property distinguishes the algorithm introduced from all other singularity-free 
algorithms known, since no artificial or constraint normalizations as well as no recursive 
or iterative procedures are necessary to conserve the rigidness of molecules. 
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C. Thermostatted dynamics 



Since the velocities appear explicitly in our approach, it is possible to introduce 
various thermostats [37-43] to simulate the canonical ensemble. Usually [26, 44], ther- 
mostatted versions allow to perform simulations with greater step sizes than those used 
within the energy-conserving dynamics. In canonical MD simulations the time evolu- 
tion should be determined in such a way to keep a fixed temperature T — (T) of the 
system, where T = 2r/ {INks) and F are the instantaneous temperature and kinetic 
energy, respectively, /cb is the Boltzmann's constant, and ( ) denotes the statistical 
averaging. The number I of degrees of freedom per particle is equal to six for arbitrary 
rigid polyatomics, except linear molecules when / = 5. 

We shall consider here a Nose- Hoover thermostat [39-41], although the extension 
to others thermostats [37, 38, 43] can also be reahzed. According to the Nose-Hoover 
technique, the thermostatted dynamics is obtained introducing the generalized fric- 
tion forces —XmVi = —Xdr/dvi. These virtual forces are added to the real ones 
and, as a result, the equations of motion for translational velocities transform into 
mdvi/dt — ii — XmVi. In our case, when orientational degrees of freedom are present 
additionally, the kinetic energy consists of both translational and rotational parts, 
^ — ^ Z^ili["^Vi^ -I- Y^a'^'^ Ja^lt^]- This requires the introduction of friction torques, 
—Xdr/dVl^^ = —XJa^a, which should be taken into account in Euler equations (4). 
The time-dependent friction coefficient varies in time with the total excess of kinetic 
energy to its canonical mean value, 

dA/di = (r-r)/(rT'), (13) 

and governs by a characteristic thermostat relaxation time r. 

In the presence of friction forces, the translational velocities can be integrated 
applying a Toxvaerd leapfrog algorithm [45]. It is based on the estimation 

Viit) = |[v,(i - I) + Vi(t + I)] + 0(h') (14) 

of on-step velocities, commonly used to calculate a kinetic part of the total energy 
at time t in microcanonical simulations and to verify the energy conservation. Then 
adding the corresponding friction term X{t)mVi{t) to the right-hand side of the first fine 
of Eq. (5) and solving the obtained equation with respect to new mid-step translational 
velocities, one finds 
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+ 1) = {[1 - imMt - 1) + hm/m}/[i + ix{t)] . (15) 

Quite analogously, in the case of rotational motion we use the estimation 

n^it) = lin^it - 1) + n,{t + 1)] + oi^) , (i6) 

for on-step angular velocities, needed to computer the friction torques —X{t)Ja^a{t). 
These torques we add to the right-hand side of Eq. (6) which now can be presented in 
the following form 

Kit + I) = {[1 - |A(t)]i^^(t - I) + /i/C^t)/ JJ/[1 + |A(t)] , (17) 

where 

Kit) = K^it) + \iJp- J,) [n^.it - Dni^it -!|) + ^Y^{t + l)ni^{t + 1)] . (is) 

The vector equation (17) constitutes a system of three equations per molecule for the 
same number of unknowns Kit + |) (note that X{t) is known from the previous time 
step). As in the case of microcanonical evaluation (A(t) = 0), the system can be 
solved iteratively, replacing initially fli{t + ^) by r2j(t — |) in all nonlinear terms which 
are collected in the right-hand side of (17). The calculated values of fli{t + |) are 
then considered as initial guess for the next iteration. The convergence of iterations is 
justified by the smallness of h which always is meet in actual MD simulations. 

Using the already defined translational and angular velocities, we evaluate the in- 
stantaneous mid-step temperature T{t+^) = T({vi(i-|- 1), fliit+ 1)}). Then equation 
(13) is integrated as follows [45]: 

X{t + h) = X{t) + h[T{t + I) - T]/{Tt^) . (19) 

The coordinates Si{t + h) and ri{t + h) are updated according to the same transforma- 
tions (see Eq. (11) and the second line of Eq. (5)) as for energy-conserving dynamics. 



D. Avoidance of iterative procedures 

An essential advantage of the approach presented lies in the fact that solutions to 
Eq. (17) (or (6)) can be obtained without applying any iterative procedures. First 
of all it is necessary to point out that the equations are nonlinear only when all the 
principal moments of inertia are different. Let us consider now this more difficult case 
(specific examples are described in Subsect. 11 F) and assume for definiteness that 
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Jx < Jy < Jz- Then the first two unknowns + |) and Qy(i + |) are the most 
fast variables which should be excluded from the iterations to increase the convergence. 

Such an excluding indeed can be carried out solving the first two (a = X, Y) 
equations of system (17) with respect to + t) and Vtyit + |). The result is 

where ^« = u_Q^,{t - |)/i/+ + [K;(t)/( 7,//+) + p^%{t - |)^(t - |)]/i, z/± = 1 ± |A(t), 
Pa = (7a/i^+, (7a = (-/^g " J-y) / C^Ja) and < //^ = -pxPv < 1/(4^'+). The last 
inequalities follow from the requirements Jq > and Jq, < J^ + Jj imposed on principal 
moments of inertia. It is worth remarking that putting formally X{t) = 0, i.e. 1/^ = 1, 
wc shall come to solutions of Eq. (6) corresponding to the microcanonical ensemble. In 
view of (20), only the third equation (a = Z) of system (17) really needs to be iterated 
with respect to one variable ^^(^ + I)- Since this variable is the most slow quantity, a 
well convergence will be guaranteed even for not so well normally behaved case as an 
almost linear body, when Jx <^ Jy < Jz- 

Finally, we shall show how to obviate the iterative solutions at all. Substituting 
result (20) into the third equation of system (17) and presenting the Z-th component of 
angular velocity in the form Q^l^ + |) = so + lead to the following algebraic equation 

ao + ai5 + 02^^ + a^5^ + 04^^ + a^5^ = (21) 

with the coefficients 

ao = (so - ez)dl - hpz[ex0Y^- + KpyOl + pxe^Y)so] , 

ai = - h^iipyOl + pxO^pz - //'so[(5so - 4^z)^+ + 2h9x9YPz]} , 

02 = /iV'[6so - 20Z + hpzOxOY + h^ii^slilOso - 6dz)] , (22) 

03 = 2/iV'[l + h^p^So{5so - 2ez)] , 

04 = h'^p^{5so — Oz) , 05 = h^jjf' , 

where '&± — 1± h?p^SQ. Equation (21) is fifth order and the corresponding solutions 
for +2) independent on parameter sq, provided the unknown 5 is precisely 
determined. However, as is well known, only algebraic equations of fourth or less orders 
allow to be solved in quadratures. 

13 



To overcome this difficulty, it is necessary to choose the parameter sq as a good 
prediction for ^^^(i + l) to be entitled to ignore the highest-order terms in the left-hand 
side of Eq. (21). The simplest choice for this can be found assuming that the nonlinear 
velocity term in the right-hand side of Eq. (18) at a = Z is timc-indcpcndcnt during 
the interval - |, t + 1] , i.e., letting (i + |)OV(i + 1) = (t - |)OV(i - |) + C>(/i) . 
As a result, one obtains 

So^ez + hpz^'xit - ^)^Y{t - I) (23) 

that represents the original values of r2^(t + |) with the second-order truncation error, 
so that 5 — 0{h?). It is easy to see that in this case the two last terms a^S'^ and a^5^ in 
the left-hand side of Eq. (21) behaves as 0{h}'^) and 0{h}'^), respectively. Taking into 
account the smallness of /i, such terms can merely be omitted without any loss of the 
precision, because they involve uncertainties of order 0{h^'^) into the desired solution 
and appear to be too small with respect to third-order truncation errors 0{h^) involved 
initially in angular velocities by the algorithm. 

Eq. (21) is now transformed into the third-order algebraic equation 

ao + ai5 + a25^ + a^S"^ = 0{h^^) (24) 

which can easily be solved analytically, 

5^ = -\a2/a^ + c-h/c + 0{h^^), (25) 
(^2,3 = -|a2/a3 - |[c - 6/c ± iv^(c + h/c)] + 0{h}'') , 

where 6= |(3aia3— a|)/a|, c = {d-\-\/W^rd?Y^^ ajudd — ^{9axa2ai—27aQal—2a\)/a^. 
Among three solutions (25), only the first one 5-i is real and satisfies the physical 
boundary condition 5\ h? when h goes to zero (the other two solutions ^2,3 are 
purely imaginary at /i — and they tend to infinity as ~ ±i//i). Thus, the desired 
Z-th component of the angular velocity is 

^'z{t+^)=so + h. (26) 

The rest two components fi^(t + |) and r2y(t + |) are reproduced on the basis of 
equalities (20). 
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E. Symplectic properties 



It can be checked readily that past and future values of all the integrated quantities 
enter symmetrically into Eqs. (5)-(7), (10), (14)-(17) and (19). Therefore, the algo- 
rithm derived is time reversible with respect to translational and rotational motions and 
within both microcanonical and canonical ensembles. In order to investigate symplec- 
tic properties, it is necessary to choose arbitrary canonically conjugated coordinates, 
express them in terms of rj, v^, Aj and fij, and, then look whether the corresponding 
volume in phase space is conserved. We accept the positions r" = rj -|- AfAa and 
momenta p" = ma{vi + Af [fij X A"]) of atomic sites to be such canonical coordinates. 

Since the translational and rotational variables are not coupled explicitly during 
our integration (in the microcanonical ensemble), the transitions from old to new 
values of the canonical coordinates, caused by varying these variables, can be con- 
sidered separately at each time step. As is well known, the translational leapfrog 
algorithm (5) is symplectic, so that the corresponding transition will be performed 
with the unit Jacobian. In its turn, the time evolution of orientational variables can 
be split into two consequent transformations. During the first one, the principal angu- 
lar velocities f2j are changed provided the orientational positions Aj remain constant, 
whereas the second transformation will correspond to the change of Aj at fixed val- 
ues of fli (a similar splitting is often used [46] to prove the symplecticity of usual 
schemes). As far as the orientational matrices A.; appear to be always orthonormal 
(detAj = 1), the effect of their changes in time is reduced simply to a rotation of 
vectors rf — rj and pf — maVj in three-dimensional space. Thus, the second trans- 
formation is evidently volume preserving. Finally, since angular velocities fij enter 
linearly into momenta p", it can be shown that the first transformation will con- 
serve the volume too provided the Jacobian jTfj. = detOj is equal to unity, where 

0, = d{n^{t + 1), nut + 1), + Dy/di^Ut - 1), ^Vit - 1), - 1)}- 

Partially differentiating each equation {a — X, Y, Z) of system (6) consequently 
over ^xit — |), f2y(t — |) and fi^(t — |), and solving the obtained three systems of 
linear equations with respect to nine elements of 0j, one finds 

_ [1 - h^jaYcrz^Sc^ + ctxCTz^y^ + crxcrY^z^} + 2h^axcrYcrznx^Y^z]t-h/2 ^^7) 
[1 - h^laYCTz^x^ + ctxctz^y'^ + (Txcty^z^} - 2h^axcrYcrznx^Y^z\t+h/2 
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Therefore, unless Jx — Jy — Jz-, the single discrete step within our algorithm is not 
volume conserving, i.e., Jciiit) — l + However, since principal angular velocities 

are Gaussian distributed in equilibrium with {Vt^VtyVt^^) = (f2^)(f2y)(r2*2-) = 0, the 
overall Jacobian Jj^ = Jl^ J^s^Xt + will fluctuate around unity. Moreover, the 
probability of deviations of from unity on some value u will decrease rapidly (as 
l/exp(~ u/h^)) with increasing u and will not depend on the total time step H 
performed. In other words, such deviations are not accumulated during the integration. 

Prom the afore said, it becomes clear that to reduce the fluctuations of JnX^) 
to zero, it is necessary to take into account the contributions (Jg — J^)^^^^.]^ = 
[W(rij) Jr2j]a of free-motion torques into the angular-velocity dynamics more precisely 
(here J denotes the diagonal matrix of principal moments of inertia). To do this, let 
us write the solutions fli{t + |) to Eq. (6) in the following formal form S{h)fli{t — |), 
where S{h) is a velocity-displacement operator. It is obvious by construction that 
]im,^^[S{h/s)Y = e^^, where L = J-^[Ki + W{Cli)3a^d/dni is a part of the Liou- 
ville operator of the system, so that dQi/dt — Lfij. Of course, the evolution operator 
e^'* does not lead, in general, to exact solutions, since the torque caused by potential 
forces is assumed to be constant, Kj = Kj(t), over the time interval — |, t + |]- Thus 
like <S{h), the propagator e^'^ generates new angular velocities with the same 0{h^) or- 
der local errors. However, the operator L allows to be decomposed into a sum -\- L2, 
and then the (Trotter) formula e = e ^ 2e^i"e ^ 2 -\-0(ji^'^ can be used to approximate 
the full propagator. Although different decompositions are possible, it is essential for 
our purposes, to decompose L in such a way to solve the problem analj^ically. 

The last requirement leads to the following decomposition 

Li = J-iK,— -, L2 = J-^W(f2,) Ji^.^ • (28) 
as li as li 

Then e^i'^fij — Q.i -\- /iJ^^Kj, whereas the propagator e^^'^/^ corresponds to a free 
rotational dynamics. As is well known, Euler equations (4) are integrated analytically 
in this case. Let \ Y.^^'^ J^^a =n&nd\ Y.^^'^ Ja^a^ = be the kinetic energy 
and square angular momenta of a body, associated with angular velocity Q, ja — 
|(-(7^(7^)-V2, e = [{Jz - Jy){M^ - mJx)lJxJYJzY''' and k = [(Jy - Jx){mJz - 
M'^)/{{Jz — Jy){M'^ — 2'HJx))Y^'^- Then, the result for transformed angular velocities 
f2' = e^^'^/^fi can be presented in the form [47]: 
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n'y = sign(nz)fi:jy£sn(£| + v) , (29) 
n'z = sign{nz)jz£dn{e^ + v) , 

where sign(a;) denotes sign of x, dn(x) = [1 — sin^(am(a;))]^/^, sn(x) = sin(am(a;)) 
and cn(x) = cos(am(a;)) are Jacobian elliptic functions [48], so that S(am(a;), = 
jam(a;) — sin^ -0] = X with S being a Legendre eUiptic integral of the first 

kind (the program code for these functions can be found in [49]). The constant v 
is determined from the initial boundary condition lim/j^o — which leads to 
V = iJsign(f2yf2^sn(u)), where v = S(arccos(f2x/(jx't^)), '^)- Note that analytical 
solutions (29) are applicable when J\4^ > 2?iJy. Otherwise, the indexes X and Z 
should be replaced between themselves in Eq. (29) and in expressions for e, k and v. 

It can be checked directly that the Jacobian of transformation f2' = e^^'^/^fi is equal 
to unity This valid also for the operator e^^'*, because a simple shift does not change 
the volume. Therefore, if the angular velocities are integrated, instead of (6), as 

12.(t + I) = e-^^te-^i V^tn.(t - I) , (30) 
our algorithm will be symplectic in rigorous meaning. 

F. Integration in specific cases 

Although the algorithm described can be applied to arbitrary rigid polyatomics, 
some simplifications are possible using special properties of the molecule. The simplest 
case is molecules with a spherical distribution of mass, when Jx — Jy — Jz = J- 
Then, no free-motion torques appear, and it is more convenient to work within the 
body- vector representation and to rewrite equations (1) and (4) in terms of angular 
velocities uJi — A^Oj in the laboratory frame, i.e., dAj/di — AiW{u)i) and JduJi/dt — 
kj — XJuji- The leapfrog trajectories for these equations are obvious: a;j(i -|- |) = 
[v_{t)uJi{t - I) + )\<ii{t)\lv+{t) and Ai{t + h) ^ Ai(t)ex.p[ipiW{<jJi)/u;il^, where 
(fi = arcsin[/iu;j/(l + ^ujf)]t+h. 

For some particular models, the orientational part of intermolecular potentials 
can be expressed using only unit vectors passing through the centers of mass of 
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molecules. The examples are point dipole interactions, when Pj = Ui/vi with Vi be- 
ing the dipole moment, or when all force sites of the molecule are aligned along Pj. 
If then additionally the condition Jx,y,z = J is satisfied (for the last example this 
can be possible when forceless mass sites are placed in such a way to ensure this 
condition), it is no longer necessary to deal with orientational matrices or quater- 
nions. In this case the equation for Pj looks as dp^/dt — W+(a;j)pj with the solution 
p.{t + h)^ exp[-^piW{uJi)/uJi]^^hPi{t). 

For molecules with the cylindric symmetry of mass distribution, the numerical 
trajectory can also be determined in a simpler manner. Let us assume for definiteness 
that Jx — Jy ^ Jz ^ 0- Then arbitrary two perpendicular between themselves axes, 
lying in the plane perpendicular to Z-th principal axis, can be chosen initially as X- 
and y-th principal orths. The corresponding solution to Eq. (17) a,t a — Z is found 
now exactly, namely, il'zit + |) = [iy_{t)n^{t - |) + j^K'z{t)]/iy+{t) (the X- and Y-th 
components are obtained automatically in view of Eq. (20)), whereas the orientational 
matrices or quaternions are computed via Eq. (11). 

A special attention should be paid on linear molecules when Jx — Jy — J ^ Jz — 
0. Each such molecule has two orientational degrees of freedom and to reproduce a 
correct dynamics by Eulcr equations it is necessary to putt formally fi^ = to exclude 
nonexisting torques caused by irrelevant rotations of the molecule around Z-axis. Then, 
it follows from Eq. (20) thatn^^(i+|) = [i/_(i)0^_y(i-|) + Jis:j^_y(i)]/i/+(i). Planar 
molecules do not present a specific case within our approach and they are handled as 
tree-dimensional bodies. 

III. MD TESTS. COMPARISON WITH PREVIOUS METHODS 

The system chosen was the TIP4P model of water [50] at a density of mN/V= 1 
g/cw? and temperature of 7=298 K. Such a system should provide a very severe test 
for rotational algorithms because of the low moments of inertia of the molecule and 
the large torques due to the site-site interactions. To reduce cut-off effects we used 
a cubic sample of iV = 256 molecules and the reaction field geometry [51]. Our MD 
programs were implemented using Fortran language and double precision throughout. 
They were executed on a Pentium-S 120 MHz personal computer at around 0.8 s per 
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time step. All runs were started from an identical well equilibrated configuration. 

We have made comparative tests on the basis of our advanced leapfrog algorithm, 
the implicit leapfrog algorithm of Fincham [26], pseudo-site formalism [23], angular- 
velocity Verlet method within matrix- and quaternion-constraint schemes [27, 28] , and 
the atomic-constraint technique [17, 18]. The results obtained for the total energy 
fluctuations £ = [{{E - {E)Y)Y''^ /\{E)\ as functions of the length {M = t/h) of the 
microcanonical (NVE) simulations are shown in Fig. 1 at four fixed step sizes, h — 
1, 2, 3 and 4 fs (a step size of 2 fs is normally used [52, 53] to simulate water within 
the atomic-constraint approach). As can be seen, the Fincham's leapfrog (marked 
simply as "leapfrog" in Fig. 1) and pseudo-site schemes appear to be unstable already 
for the least time step considered and they lead to the worst energy conservation. 
Much more stable trajectories are produced by the velocity- Verlet integrator within 
quaternion- and matrix-constraint schemes which exhibit similar equivalence in the 
energy conservation. But the results are rather poor at moderate and great time steps 
(/i > 3 fs). Only the atomic-constraint and our advanced leapfrog algorithms can be 
related to long-term stable schemes. 

We mention that the computation of total energy E — F+U (consisting of kinetic F 
and potential U parts) at time t within the leapfrog framework requires the knowledge 
of on-step velocities. These velocities can be calculated using usual estimators (14) and 
(16). The corresponding curves of dependencies ^(jV) are labeled by "1" in Fig. 1. They 
are almost identical to those obtained within the atomic-constraint technique. Note 
that 0{h'^) uncertainties, arising in evaluations (14) and (16), are in the self-consistency 
with second-order global errors appearing during the leapfrog integration. Despite this, 
an additional portion is involved into the main 0{h'^) term of accumulated errors for 

increasing the total energy fluctuations with no relation to the real accuracy of the 
computed trajectory [54]. Although various more accurate estimators are available 
[55], we have established that the following four-point symmetric scheme 

Y{t) = ^ [ - V(i - f ) + 9V(i - I) + 9V(i + I) - V(i + f )] + 0{h^) , (31) 

where V = {vj, f2j}, leads to the best energy conservation within our leapfrog algo- 
rithm. It is understood, of course, that mid-step velocities entering into the right-hand 
site of Eq. (31) are already defined quantities. Thus, the computation of the total 
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energy E at time t becomes possible, when the velocity step t + ^ has been passed 
and the potential energy U (t) (known at this stage already at t + h) has been taken 
from memory. The corresponding dependencies of S on Af are plotted in Fig. 1 by the 
boldest (lowest lying) curves marked as "2". In this case, the total-energy fluctuations 
decrease about in 1.5 times with respect to the usual two-point scheme. 

The quaternion and principal-axis representations of the advanced leapfrog algo- 
rithm conserved the energy approximately with the same accuracy. For this reason, the 
result concerning principal- axis variables is not plotted in Fig. 1 to simplify the graph. 
The leapfrog trajectories were generated applying quasianalytical solutions (Eqs. (20) 
and (26)) for angular velocities. The solutions obtained by means of iterations of 
Eq. (6) were calculated also for comparison. No deviations between the both results 
have been identified up to /i = 6 fs. They differed on each step by uncertainties of 
order round-off errors only, so that the free-of-iteration scheme appears to be in an 
excellent accord. No shift of the total energy and temperature was observed during 
the advanced leapfrog trajectories at /i < 5 fs over a length of 10 000 time steps. The 
deviations from unity of the overall Jacobian Jj^ (see Eq. (27)) never exceed about 5% 
(at < 4 fs). 

Finally, the rigorously symplectic version (30) has also been examined (the longest- 
dashed curves in Fig. 1) within the four-point scheme (31). As we can see, this version 
does not lead to improvements in energy conservation, despite the fact that the free- 
motion propagator e^^^^'^ is evaluated exactly. This is so because for water at the given 
thermodynamic point, the free- motion torques are much smaller in amplitude with 
respect to the torques caused by interactions. An increased conservation of energy can 
be expected for systems at high temperatures, where the free-motion contributions into 
the torques become dominate. Otherwise, the analytical version (29) is not generally 
recommended because it requires the calculation of somewhat time-consuming elliptic 
functions (although we did not observe any considerable decreasing time efficient, given 
that near 95% of the total computer time were spent to evaluate pair interactions). 
In the limiting case when the potential-force torques are absent at all, symplectic 
solutions (29) will lead to exact results with automatic preservation of energy and 
angular momenta. The symplectic version should also be used in situations where a 
precise conservation of the volume in phase space at each time step is very important. 
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It is worth remarking that the variant Li <-> L2 of decomposition (28) in transformation 
(30) is also acceptable in view of the symplecticity. We have estabhshed, however, that 
such a decomposition leads to worse energy conservation in our simulations (probably 
because of (K^^) > ([W(rii) jrii]^)). 

To properly reproduce features of NVE ensembles, it is necessary for the ratio T of 
£ to the fluctuations U of potential energy to be no more than a few per cent. We have 
obtained the following levels of £ (within the two-point scheme) at the end of these 
trajectories: 0.0016, 0.0065, 0.015, 0.029, 0.049 and 0.10 %, corresponding to T « 0.29, 
1.2, 2.7, 5.2, 8.7 and 18 % at h = 1, 2, 3, 4, 5 and 6 fs, respectively (U ^ 0.56% for the 
system under consideration). Therefore, a step size of 4 fs is still suitable for precise 
calculations. The greatest time steps 5 fs and 6 fs can sometimes be acceptable when 
the precision is not so important, for example, for the equilibration of configurations. 
The ratio T in the interval h < 5 fs can be fitted with a great accuracy to the function 
Ch^ with the coefficient C ~ 0.29 % fs~^. This is completely in line with the square 
growth of global errors appearing during the integration by Verlet-type integrators [54] . 

To verify the angular-velocity approach in more detail, in our NVE and canonical 
(NVT) ensemble simulations we measured besides the total energy and temperature 
some other relevant functions of the system, namely, specific heat at constant vol- 
ume, mean-square forces and torques, oxygen-oxygen and hydrogen-hydrogen radial 
distribution functions (RDFs). Center-of-mass (CM) and angular- velocity (AV) time 
autocorrelation functions (TAFs) were also found. Orientational relaxation was studied 
by evaluating the molecular dipole-axis (DA) autocorrelations. The NVE simulations, 
performed within the atomic-constraint technique at h = 2 fs, was considered as a 
benchmark against which other algorithms and step sizes are to be compared. First 
of all, to finish the discussion with the NVE integrators, we report that deviations in 
all the measured functions with respect to their benchmark values were in a complete 
agreement with the corresponding relative fluctuations T. For example, the results 
obtained with the help of the advanced leapfrog algorithm at /i = 2 fs were indistin- 
guishable from the benchmark ones. At the same time, they differed as large as around 
5%, 10% or even 20% with increasing the time step to 4 fs, 5 fs or 6 fs, respectively, 
however, these differences were smaller than for other integrators. Therefore, there is 
a little point in pursuing the energy-conserving algorithms to time steps larger than 4 
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fs because the deviations become evident. 

In the case of NVT simulations, the investigated quantities were less sensitive to 
the step size increasing. These simulations have been performed using the Fincham 
implicit algorithm [26] (within the Brown and Clarke thcrmobath [37]) as well as the 
advanced leapfrog integrator within the Nose-Hoover thermostat (Sec. II C) applying 
the free-of- iteration scheme (Eqs. (20) and (26)). The thermostat relaxation time r was 
chosen to be 1 ps, (i.e., h <^ r <^ Afh) and the friction coefficient X{t) was putted to be 
zero at the very beginning (t = 0) of the simulations. The RDFs and CM, AV and DA 
TAFs calculated during the NVT integrations at different step sizes {h = 2-10 fs) are 
shown in Fig. 2 in comparison with the benchmark results. These functions at h — 4 
and 6 fs coincided completely with those corresponding to /i = 2 fs and they are not 
shown in subsets (a)-(b) to simphfy the presentation. As can be seen from Fig. 2, the 
RDFs remain practically the same at < 8 fs in the case of the advanced integrator. 
The deviations of all the correlation functions from the benchmark within the usual 
implicit algorithm are clearly larger. The orientational correlation function (see subset 
(d)) is to show a systematic discrepancy in this case. For instance, these deviations at 
h — A is are as big as those obtained during the advanced leapfrog integration ai h — 
8 fs. Therefore, the last approach allows a step size approximately twice larger than 
the usual implicit integrator. We conclude, therefore, that the thermostatted advanced 
leapfrog integrator allows to be used up to a time step of 6 fs, given that then there is 
no difference in RDFs, while CM, AV and DA TAFs are also close. 



IV. CONCLUDING REMARKS 



We have formulated a new approach for numerical integration of the equations 
of motion for systems with interacting rigid bodies. Unlike other standard meth- 
ods, the principal angular velocities are involved directly into the integration within 
this approach. The algorithm derived is categorized as a rotational leapfrog, since 
the variables saved are mid-step angular velocities and on-step orientational positions. 
The orientations can be expressed in terms of either quaternions or entire rotational 
matrices. The interpolation of velocity- and orientation-dependent quantities to the 
corresponding middle time points was carried out using a simple averaging over the 
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two nearest neighboring values that is in the spirit of the leapfrog idea as well. As 
a result, the following significant benefits have been achieved: (i) the exact conserva- 
tion of rigid structures appears to be an intrinsic feature of the algorithm, and (ii) all 
the evaluations arc performed analytically in both NVE and NVT ensembles without 
involving any iterative procedures. 

It has been shown on the basis of an actual computer experiment on water that the 
algorithm presented exhibits better energy conserving properties than those observed 
in all other rigid-motion integrators known. The algorithm can easily be implemented 
for arbitrary rigid bodies and substituted into existing program codes. It seems to 
be a good alternative to the atomic-constraint method in the case of long term MD 
simulations of systems with rigid polyatomic molecules. 

The approach introduced can be developed to perform a MTS integration within the 
RESP technique. This question will be considered in our further studying. Some ideas 
of the MTS integration have already be used in this paper to derive symplectic versions 
of the advanced leapfrog algorithm. The reformulation of the RESP approach in gen- 
eralized coordinates containing translational and orientational variables explicitly may 
lead to significant simplifications when selecting efficient reference system propagators. 
For example, translational motions (which are much slower than rotational ones for the 
most of liquids) can be decomposed directly within such coordinates. The most noto- 
rious demonstration is very fast rotations of rigid bodies in the inertial-motion regime. 
Then for molecules with a spherically symmetric distribution of mass we merely ob- 
tain Vj — const and fij = const within the molecular approach. To reproduce such a 
behavior within the atomic technique, it is necessary to integrate (with a very small 
step size) the equations of motion in the presence of rapidly varying strong constraint 
forces. The decomposition of rotational propagators containing interactions into slow 
and fast parts can also be done easily splitting in an appropriate way the free-motion 
and potential-force torques. 
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Figure captions 



Fig. 1. The total energy fluctuations as functions of the length of the NVE simu- 
lations on the T1P4P water, performed in various techniques at four fixed time steps: 
(a) 1 fs, (b) 2 fs, (c) 3 fs and (d) 4 fs. 

Fig. 2. Oxygen-oxygen (0-0) and hydrogen- hydrogen (H-H) radial distribution 
functions (a), center-of-mass (b) and angular- velocity (c) time autocorrelation func- 
tions, and orientational relaxation (d), obtained in the NVT simulations of the TIP4P 
water. The results corresponding to the step sizes h — 2,S and 10 fs are plotted by bold 
solid, short-dashed and thin solid curves, respectively. Additional long-short dashed 
and dashed curves in (c)-(d) correspond to the cases of h = A and 6 fs. The sets of 
curves related to the usual implicit and advanced leapfrog algorithms are labeled by 
"U" and "A", respectively (the result of the usual implicit algorithm is not included in 
(c) to simplify the graph). The benchmark data are shown as open circles. Note that 
the advanced-algorithm curves are indistinguishable in (d) at /i = 2, 4 and 6 fs. 
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